/*
  Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

  This file is part of Sunrise.

  Sunrise is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 3 of the License, or
  (at your option) any later version.

  Sunrise is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.
  
*/

// $Id$

#include<iostream>
#include "CCfits/CCfits"
#include "counter.h"
#include "polychromatic_grid.h"
#include "optical.h"
#include "scatterer.h"
#include "constants.h"
#include "dust_model.h"
#include "emission.h"
#include "emission-fits.h"
#include "full_sed_emergence.h"
#include "full_sed_grid.h"
#include "emergence-fits.h"
#include "mcrx-units.h"
#include "grid_factory.h"
#include "grid-fits.h"
#include "shoot.h"
#include "fits-utilities.h"
#include "boost/shared_ptr.hpp"
#include "density_generator.h"
#include "preferences.h"
#include "model.h"
#include "wd01_Brent_PAH_grain_model.h"
#include "wd01_grain_model.h"
#include "xfer.h"
#include "ir_grid.h"
#include "equilibrium.h"

using namespace blitz;
using namespace mcrx;
using namespace CCfits;
using namespace std; 

typedef local_random T_rng_policy;
typedef scatterer<polychromatic_scatterer_policy, T_rng_policy> T_scatterer;
typedef T_scatterer::T_lambda T_lambda;
typedef T_scatterer::T_biaser T_biaser;
typedef dust_model<T_scatterer, cumulative_sampling, T_rng_policy> T_dust_model;
typedef full_sed_emergence T_emergence;

typedef full_sed_grid<adaptive_grid> T_grid; 
typedef ir_grid<adaptive_grid, cumulative_sampling, local_random> T_ir_grid;
typedef emission<polychromatic_policy, T_rng_policy> T_emission;


class scattered_intensity_factory: public grid_factory<T_grid::T_data> {
public:
  const T_densities rho_;
  const T_float scatterer_z_;
  const int n_lambda;
  scattered_intensity_factory (const T_densities& rho, T_float z, int nl):
    rho_(rho), scatterer_z_ (z), n_lambda(nl) {};
  
  virtual bool refine_cell_p (const T_cell& c, int level) {
    return level <0;};
  virtual bool unify_cell_p (T_grid::const_local_iterator begin,
			     T_grid::const_local_iterator end) {
    return false;};
  virtual T_data get_data (const T_cell& c) {
    return T_data (rho_, vec3d(0,0,0), T_lambda(n_lambda));
  }
    
  virtual int n_threads () {return 1;};
  virtual int work_chunk_levels () {return 5;};
  virtual int estimated_levels_needed () {return 0;};
};


void scattered_intensity (T_float rho=0.01, int grid_n=1,
			  int n_forced=1, int seed = 0)
{  
  //const T_float opacity = .01;
  const T_float emitter_z = -50;
  //  const T_float scatterer_z = 50;
  //const T_float emitter_z = -interaction_d;
  const T_float scatterer_z =0;
  const T_float grid_extent = 1;
  //const int grid_n = 1;
  const T_float grid_cell_size = 2*grid_extent/grid_n;
  const T_float cell_volume = pow (grid_cell_size, 3);
  const int n_cameras = 17;
  const T_float distance = 10000;
  const int n_threads=8;
  const int n_lambda=1;
  const int n_rays=10000000;

  T_lambda lambda(n_lambda);lambda=1;

  mcrx::seed(seed);
  T_unit_map units;
  units ["length"] = "kpc";
  units ["mass"] = "Msun";
  units ["luminosity"] = "W";
  units ["wavelength"] = "m";
  units ["L_lambda"] = "W/m";

  // load grain model
  T_lambda kk(n_lambda); kk=1;
  T_lambda aa(n_lambda); aa=.9;
  T_lambda gg(n_lambda); gg=.5;
  T_dust_model::T_scatterer_vector sv;
  sv.push_back(boost::shared_ptr<T_scatterer> 
	       (new simple_HG_dust_grain<polychromatic_scatterer_policy, mcrx_rng_policy> (kk,gg,aa)));
  T_dust_model model (sv);
  model.set_wavelength(lambda);
  
  const T_float albedo = model.get_scatterer(0 ).albedo()(0);

  T_densities rhov(n_lambda);rhov=rho;
  scattered_intensity_factory factory (rhov, scatterer_z, n_lambda);
  // first create the adaptive grid
  boost::shared_ptr<adaptive_grid<T_grid::T_data> > tempgrid
    (new adaptive_grid<T_grid::T_data>  
     (vec3d (-grid_extent,-grid_extent,-grid_extent+ scatterer_z ), 
      vec3d (grid_extent,grid_extent,grid_extent + scatterer_z),
      vec3i (grid_n, grid_n, grid_n), factory));
  
  // and then the full_sed_grid
  T_grid gr(tempgrid, units);

  //T_grid g(structure_hdu, data_hdu, model, 500);
  cout << "setting up emission" << endl;
  auto_ptr<T_emission> em;
  //em.reset(new pointsource_emission<polychromatic_policy, T_rng_policy> 
  //	   (vec3d (0., 0., emitter_z), lambda));
  em.reset(new laser_emission<polychromatic_policy, T_rng_policy> 
	   (vec3d (0., 0., emitter_z), vec3d(0,0,1), lambda));
  
  auto_ptr<T_emergence> e(new T_emergence(n_cameras,1, distance, 3*grid_extent, 1));

  std::vector<T_rng::T_state> states= generate_random_states (n_threads);
  T_rng_policy rng;

  xfer<T_dust_model, T_grid, T_rng_policy> x (gr, *em, *e, model, rng, 1e-2, 0,n_forced,10.);
  long n_rays_actual = 0;
  dummy_terminator t;
  T_biaser b(0);
  shoot (x, scatter_shooter<T_biaser> (b), t, states, n_rays,
  	 n_rays_actual, n_threads, true, n_threads,0, "");

  T_float normalization = 1./n_rays_actual;

  //FITS output("!scattered.fits", Write);
  //e->write_images(output,normalization,units, "", false, false);

  // analyze flux in images
  T_float four_pi = 0;
  int n = 0;
  const T_float d = scatterer_z-emitter_z;
  const T_float flux_at_dust =n_rays_actual/(4*grid_extent*grid_extent);
  const T_float scattered_luminosity =
    flux_at_dust*rho*pow(2*grid_extent,3)*albedo;
  for (T_emergence::iterator i = e->begin(); i != e->end(); ++i) {
    const T_emergence::T_camera::T_image image =
      (*i)->get_image ();
    const T_float scatter = sum (image);
    if ((i != e->begin()) && (i != e->end()- 1)) { 
      four_pi += scatter;
      ++n;
    }
    const T_float theta =(*i)->get_theta();
    const T_float phase = model.get_scatterer(0 ).phase_function().distribution_function(vec3d(0,0,1),vec3d(sin(theta),0,cos(theta)))(0);

    const T_float theoretical = scattered_luminosity*phase/(distance*distance);


    cout << (*i)->get_theta()
	 << "\t" << scatter
	 << "\t" << theoretical << '\t' << scatter/theoretical << endl;

  }
  /*
  cout << "TOTAL QUANTITIES\n";
  const T_float total_scattered = four_pi/n*(4*constants::pi*distance*distance);
  cout << "\ntotal scattered luminosity " << total_scattered << endl;
  cout << "should be " << scattered_luminosity << endl;
  cout << "discrepancy is " << (total_scattered - scattered_luminosity)/scattered_luminosity << endl;
  */
}

void test_scattered_intensity ()
{
  //cout << "Testing opacity" << endl;
  //for (T_float o = 0.01;  o <= 10; o*= 10) {
  //cout << "opacity: " << o << endl;
  //  scattered_intensity (o, 50, 1);
  // }
  cout << "tau=0.001, grid_n=10, n_forced=0\n";
  scattered_intensity (0.001, 10, 0, 42);
  cout << "tau=0.001, grid_n=1, n_forced=0\n";
  scattered_intensity (0.001, 1, 0, 42);
  cout << "tau=0.001, grid_n=10, n_forced=1\n";
  scattered_intensity (0.001, 10, 1, 42);
  cout << "tau=0.001, grid_n=1, n_forced=1\n";
  scattered_intensity (0.001, 1, 1, 42);

  cout << "tau=1.0, grid_n=10, n_forced=0\n";
  scattered_intensity (1.0, 10, 0, 42);
  cout << "tau=1.0, grid_n=1, n_forced=0\n";
  scattered_intensity (1.0, 1, 0, 42);
  cout << "tau=1.0, grid_n=10, n_forced=1\n";
  scattered_intensity (1.0, 10, 1, 42);
  cout << "tau=1.0, grid_n=1, n_forced=1\n";
  scattered_intensity (1.0, 1, 1, 42);
  cout << "tau=5.0, grid_n=1, n_forced=0\n";
  scattered_intensity (5.0, 1, 0, 42);
  cout << "tau=5.0, grid_n=1, n_forced=1\n";
  scattered_intensity (5.0, 1, 1, 42);
}

int main (int, char**)
{
  test_scattered_intensity();
}
